Propagation of fluctuations in interaction networks governed by the law of mass action. 
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Using an example of physical interactions between proteins, we study how perturbations propagate in inter- 
connected networks whose equilibrium state is governed by the law of mass action. We introduce a comprehen- 
sive matrix formalism which predicts the response of this equilibrium to small changes in total concentrations 
of individual molecules, and explain it using a heuristic analogy to a current flow in a network of resistors. Our 
main conclusion is that on average changes in free concentrations exponentially decay with the distance from 
the source of perturbation. We then study how this decay is influenced by such factors as the topology of a 
network, binding strength, and correlations between concentrations of neighboring nodes. An exact analytic 
expression for the decay constant is obtained for the case of uniform interactions on the Bethe lattice. Our 
general findings are illustrated using a real biological network of protein-protein interactions in baker's yeast 
with experimentally determined protein concentrations. 
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Equilibria in a broad class of microscopically reversible 
processes whose direct and reverse rates are proportional to 
the product of concentrations of participating substances is 
described by the Law of Mass Action (LMA). It has been rig- 
orously proven that the unique steady state of such processes 
is completely defined by the set of initial concentrations and 
reaction constants Recently, it has become popular to vi- 
sualize large sets of interacting substances as networks, where 
nodes and links correspond to the reactants and their propen- 
sity for pairwise interactions. One of the best-studied exam- 
ples of such networks is that formed by all protein-protein 
physical interactions (pairwise bindings) in a given organ- 
ism. The LMA determines the equilibrium free concentrations 
of individual proteins as well as those of hetero- and homo- 
dimers. A surprising feature observed in virtually all recent 
large-scale studies of these networks in a wide-ranging vari- 
ety of biological organisms is their globally connected topol- 
ogy. Indeed, most pairs of protein nodes are linked to each 
other by relatively short chains of interactions. Total (bound 
plus unbound) concentrations of individual proteins are sub- 
ject to both stochastic fluctuations due to the noise in their 
production and degradation as well as to systematic changes 
in response to external and internal stimuli. Such localized 
perturbation changes bound and free concentrations of imme- 
diate neighbors of the perturbed protein which in their turn 
influence their neighbors, etc. Thus a fluctuation in the to- 
tal concentration of just one reactant to some degree affects 
free concentrations of all substances in the same connected 
component of the network. The propagation of fluctuations 
far away from their source presents a great threat of undesir- 
able cross-talk between different functional processes, simul- 
taneously taking place in an organism. Thus it is important 
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to understand whether and how this propagation gets attenu- 
ated and under what conditions it is minimized. There also 
exist several anecdotal cases when changes in free concentra- 
tions propagating beyond the immediate neighbors of a per- 
turbed protein are used for a meaningful biological regula- 
tion/signaling. Then a relevant question is under what condi- 
tions the propagation of the signal in a desirable direction is 
the least attenuated. 

In this work we study how localized perturbations such 
as changes of concentrations of individual proteins affect the 
binding equilibrium at all nodes of a protein interaction net- 
work. Our main conclusion is that under a broad range of con- 
ditions such perturbations exponentially decay with the net- 
work distance away from the perturbed node. Luckily, this 
makes protein binding networks poor conduits for indiscrimi- 
nately propagating fluctuations which would have led to unde- 
sirable cross-talk between biologically distinct pathways. On 
the other hand, under carefully selected conditions, a pertur- 
bation can propagate relatively far with a minimal attenuation. 
We first consider a general case of propagation of small per- 
turbations in a network of an arbitrary topology, concentra- 
tions, and dissociation constants. Several simpler analytically 
tractable case studies follow. 

Consider a network of pairwise interactions between N dis- 
tinct types of proteins (or any other molecules for that matter), 
where each protein corresponds to a vertex. The existence of a 
link between vertices i and j means that substances i and j in- 
teract to form a bound complex (hetero- or homo- dimer) ij. 
Throughout this paper we would consider only such dimers 
and ignore the existence of multi-protein complexes consist- 
ing of 3 or more proteins. However, we checked |2] that our 
main results could be easily extended to an arbitrary compo- 
sition of complexes or even to a situation in which individual 
nodes combine with each other in reversible chemical reac- 
tions obeying the LMA. 

The LMA dictates that free concentrations of proteins Fi 
and those of dimers Dij obey 
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where fc^ is the corresponding dissociation constant. Adding 
the mass conservation one obtains the following system of 
equations that relates total concentrations C; of proteins to 
their free concentrations F, , 



Ci — Fi 



E 



F,F, 



(2) 



Here and below the notation means a sum over all ver- 

tices j that are the network neighbors of the vertex i. In a 
general case of 4 or more interconnected interacting pairs this 
system of non-linear equations allows for only numerical so- 
lution. One particularly convenient computational method in- 
volves rewriting the Eq. ^ as 



F, = 



1 + E,^. F,/h 
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and successively iterating it starting with Fi ~ Ci until the 
LMA is satisfied with a desired precision. The proof presented 
in Ref. guarantees the uniqueness of the solution found 
this way. Now consider a change in free concentrations 5Fj 
induced by a small perturbation of total concentrations 6Ci. 
Linearization of Eq. (|2]l yields 
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with matrix A ^ given by 
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Here Aij is the adjacency matrix of the network. In case of 
multi-protein complexes consisting of more than 2 different 
proteins Dij should be replaced with the total concentration 
of all complexes containing both i and j among their con- 
stituents. It follows from Eqs. ( 1415b that if the change in total 
concentration was limited to just one node 0, the induced rel- 
ative change of free concentration of any other node i 7^ 
satisfies 



5F, _ ^ SFj 
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This equation shows that changes in free concentrations on 
nearest neighbors tend to be of the opposite sign. Also, 
since Ejv^i^u/^j ~ 1 ~ Fi/Ci < 1, the absolute mag- 
nitude of perturbation \SFi/Fi\ on any node away from the 
source is less or equal than its maximal value among its neigh- 
bors: maxj^i \dFj/Fj\. Bonds with higher Z^ij/C,; are better 
transmitters of perturbations from node j to node i. Note, that 
this quantity is non-symmetric: The transmission along any 
particular edge is directional with preferred direction pointing 
from the higher total concentration to lower one. 

Inverting ^ one obtains the desired expression for the lin- 
ear response of any free concentration to an arbitrary pertur- 
bation in total concentrations: 
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FIG. 1: A) The average magnitude of normalized changes in free 
concentrations \5Fi/Fi\ per unit of 5Co/Co plotted as a function of 
the distance Uo from the perturbed node 0. The propagation of per- 
turbations (Eq. [7} was computed in a highly curated yeast PPI net- 
work 1 4] with real-life total concentrations of individual proteins |3] 
and dissociation constant Kd increasing in 10-fold increments start- 
ing from InM (or 34 molecules/cell) and up to 0.1 mM (the steepest 
decaying curve). 

B) The exponential decay constant ^(Kd) obtained from the best fit 
in the form A expl—'y {Kd)lio] to curves shown in the Panel A with 
the real-life concentrations (filled circles), and for randomly reshuf- 
fled concentrations (open diamonds). 



As an application of the general theory outlined above we 
calculated the free concentrations Fi and the linear response 
matrix A^^ in a realistic protein network of baker's yeast 
Saccharomyces cerevisiae. A highly-curated set of protein- 
protein physical interactions from the BIOGRID dataset \^ 
was further restricted to include only interactions which were 
reported in at least two publications. Total concentrations of 
proteins for yeast grown in rich growth medium conditions 
were taken from a genome-wide experimental study |5]. The 
resulting dataset consists of 4185 interactions among 1740 
proteins with total concentrations ranging between 50 and 1 
million molecules/cell with median ~ 3000 molecules/cell. In 
the absence of genome-wide information regarding the value 
of dissociation constants in our simulations we assume them 
all to be the same fcy = Kd- We also tried simulations in 
which kij are randomly drawn from a particular probability 
distribution. In all of our simulations we observed that the 
magnitude of relative changes in free concentrations exhibits a 
universally exponential decay with the network distance from 
the source of perturbation (Fig. [TJ\). That is to say, on aver- 
age, matrix elements of (A^^)^^ exponentially decrease as a 
function of distance between i and j. Further progress in ex- 
perimental methods will evidently lead to discovery of new 
protein interactions and more precise values of protein con- 
centrations and dissociation constants. However, such im- 
provements will not affect our general conclusion about the 
exponential decay of perturbations. 

While the linear response approximation (|7]) describes in- 
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finitesimally small perturbations, the response to a finite per- 
turbation could be obtained numerically by repeating itera- 
tions of ^ with perturbed total concentrations and compar- 
ing the resulting free concentrations to the wild-type ones. 
We simulated the response of the system to a 20% decrease 
(which is roughly the range of intrinsic fluctuations JBI]) as 
well as to a complete elimination (which is experimentally re- 
alizable as a gene knock-out or inactivation) for each of the 
proteins in our dataset. It was found that even in the latter 
extreme case the linear response approximation works rather 
well. The exponential decay of perturbation was found to be 
identical to that calculated using the Eq. (O and the overall 
magnitude of changes in free concentration is comparable to 
that predicted from the linear response. 

To develop a better understanding of this matrix formal- 
ism, we first consider the case when the underlying network 
is bipartite (but not necessarily acyclic). To take into account 
the natural sign- alternation of 6Fi / Fi on immediate neighbors 
in the network (see Eq. we introduce the new variables 
(j>i = {—l)^'5Fi/Fi where Xi is on one sublattice and 1 on 
the other. This allows us to rewrite the Eqs. ( |4l5l l as 

5a = Y,'^v^3 ■ (8) 

Here 6Ci = {—l)^'SCi and & is given by 

aij = -A,jD,j + D,k + F^j . (9) 

In the situation when SCi = 6ioSCo (i.e. when the pertur- 
bation is limited to a single node 0), the Eqs.( |8l9t can be in- 
terpreted as describing "electric potentials" (j>i in the network 
of resistors with resistances i?^ = 1/aij = l/Dij subject 
to the injection of the current SCq at the node 0. Each node 
is also shunted to an auxiliary "ground node" with potential 
4>G = by the resistance Ric = l/Fi. Potential gradients 
along edges 0, - 0, = {-1)''^SFJF, - SFj / Fj = 

{—l)^'6Dij/Dij determine relative (dimensionless) changes 
in concentrations of heterodimers, while currents = 
{(j>i — <j)j)/Rij = {—l)^^dDij - the absolute (dimensional) 
changes. Similarly, currents to the ground lio — (j)i/RiG = 
{—ly^'SFi are equal to changes in free concentrations of pro- 
teins. As in resistor networks, the Kirchoff Law here fol- 
lows from the mass conservation which states that everywhere 
the total current flowing out of node i - lio + ^ij = 

i-ir^iSF, + Y.^^j SD,,) = i-lY-SC, = sc., is equal to 
the external current SCi — S^SCo of changes in total concen- 
trations. 

The interpretation of the free concentrations Fi as "shunt 
conductivities" leaking the "current" to the ground means that 
the smaller they are, the weaker is the decay of both currents 
SD and SF with the distance. Since stronger binding gener- 
ally decreases free concentrations of all proteins, it naturally 
reduces the rate of decay of perturbations (visible in the Fig. 
[TJ\). However, the exponential decay constant jiKd) appears 
to saturate around 2.25 as Kd (Fig. [IB)- 
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FIG. 2: The propagation ratio /ii^s defined in the text in a linear 
chain Ci — > C2 ^ C3 with strong binding (k — InM) as a function 
of concentration of the intermediate substance C2. All curves have 
the same concentration of the target protein C3 — l/iM and different 
concentrations of the source of the signal Ci = 2/iM (dashed line), 
Ci = 1/iM (solid line), and Ci = 0.5/iM (dot-dashed line). In 
the strong binding limit the propagation ratio reaches its maximum 
when all three substances are completely bound, C2 = Ci + C3. It 
also follows that the peak propagation is better from high to low total 
concentration i.e. when Ci > C3. 



This saturation is easy to understand. Indeed, consider the 
most ideal scenario in which all free concentrations Fi are 
very small and thus the "current" SD is approximately con- 
served (loss to the ground is negligible.) The exponential 
growth in the number of neighbors Af„(Z) ~ {{d{d—l))/{d)y 
as a function of distance I from the perturbation source means 
that even in this ideal setup the average current at distance 
I would be proportional to 1/Nn{l) and thus exponentially 
small. The same rate of decay describes the "potentials" 
SFi/F, as well. 

However, it is important to emphasize that the ideal sce- 
nario outlined above is almost never realized with real-life 
concentrations. Indeed, the limit of infinitely strong binding 
kij does not make all but only some of free concentra- 
tions Fi to go to zero. One can see it clearly already for two 
interacting proteins. When their concentrations Ci and C2 are 
not equal to each other, in the strong binding limit ki2 
the free concentration of the more abundant protein (say 1) 
remains nonzero Fi ^ Ci — C2, while the free concentration 
of its less abundant partner F2 0. Consider another simple 
example of a chain of three proteins with initial concentrations 
Ci, C2 and C3 reacting to form dimers 1 — 2 and 2 — 3 with 
the same dissociation constant k. In this case one could still 
analytically calculate all free and bound concentrations. The 
logarithmic derivative /ii^3 = (9_F'3/9Ci)(Ci/-F'3) quanti- 
fies the propagation of perturbation of the node 1 through this 
3 -node channel and in the strong binding limit has a maxi- 
mum around C2 — Ci + C3 (see Fig|2]l i.e. when all three 
substances are completely bound: ^1,^2,^3 0. In a gen- 
eral case of a PPI network of arbitrary topology the only sit- 
uation in which free concentrations of all proteins would ap- 
proach zero as kij ^ is when their total concentrations Ci 
are proportional to their degrees di. For a given topology of 
the network such concentration setup has the slowest decay of 
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perturbations. 

Most real-life PPI networks are characterized by a positive 
correlation between total concentrations of interacting pro- 
teins. In the yeast network used in this study we observed this 
effect to be present and highly statistically significant: (the 
Spearman rank correlation coefficient of 0.27 with P- value of 
10^^^.) Such correlation improves the balance between total 
concentrations of interacting nodes and thus somewhat low- 
ers the average free concentration of proteins compared to a 
case where this correlation is absent. Based on this we ex- 
pect that real protein-protein networks would be more prone 
to propagating perturbations than their counterparts in which 
concentrations of proteins are randomly reshuffled and thus 
the correlation between concentrations of interacting nodes is 
destroyed. This theoretical expectation was indeed verified in 
Fig. [TJj (compare filled circles and open diamonds for the 
network with real concentrations and the reshuffled ones). 

Most real-life PPI networks are not bipartite. Fortunately, 
due to a relative sparsity of links the number of odd-length 
loops in them is small and our resistor network analogy pro- 
vides a reasonable approximation. For any starting point of 
perturbation the optimal way to define sign-alteration in 
variables (pi = {—iy^'6Fi/ Fi is by using Xi = lio (here ko is 
the distance from the source of perturbation to the node i.) 
The majority of links would connect nodes with opposite 'par- 
ities', while the remaining non-bipartite links could be treated 
as a small but important correction to the ideal case. Like 
shunt conductivities to the ground, they contribute to the dis- 
sipation of the "current". Indeed, if a link z ^ j is of this 
anomalous kind, its contribution to the current leaving nodes i 
and j are equal to each other and given by Dij {4^ + 0j ). One 
example of such anomalous (non-bipartite) links is given by 
homodimers (see |3].) In general, these anomalous links lead 
to the loss of the total current from the system and thus tend 
to dampen the propagation of perturbations. 

In what follows we analytically investigate a simple exam- 
ple of a bipartite network - the Bethe lattice - where all dis- 
sociation constants are the same, fcy = k, and each vertex 
has the same number of interaction partners (degree) di = d. 
When total concentrations of all proteins are also identical 
Ci — C, their equilibrium free concentrations are given by 



be rewritten as 



F,=F = 




(10) 



Small deviations 6Fi at distance I from the perturbation source 
at node should satisfy the Eq. Q, which in this case could 



(d + -)6Fi = (d - l)SFi+i + SFi^i 

r 



(11) 



This equation has an exponentially decaying solution 5Fi = 
SFq a', where 



d + k/F 



'd + k/F 
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As expected, — 1 < A < which means that perturbations 
sign alternate and exponentially decay as a function of I. In a 
strong binding limit the combination of Eqs. ( I12I10I ) yields 



1 / 1 dk\ ^^,dk. 



(13) 



For a linear chain of proteins (d — 2), the complete solution 
in terms of C and k looks particularly elegant: 



(l + 8C/fc)i/'^-l 
"(l + 8C/fc)i/4 + i- 



(14) 



In the limit of strong binding, a perturbation in a linear chain 
propagates indefinitely, |Aid| 1. Using this approach, 
it also turns possible to find the exact solution for the de- 
cay exponent in the linear chain (d — 2) with oscillating 
total concentrations of proteins: Ci = C[l + (— l)'a]. Re- 
sponse of the even- and odd-numbered vertices has different 
amplitudes yet decays with the same exponential coefficient 
A± - -[4xVT^ - V2/( l + 4x)-27^]/[l + 4x - /], 
where x ~ C/k and / = y^lGa^x^ + 8% + 1. Evidently, 
|A±(C/A:, a)\ < \l{C /k) with equality being achieved only 
for a — Q. Thus, as was argued in a general case, any variation 
in Ci /di results in a larger average unbound concentrations Fi 
and thus in a faster decay of perturbations. 
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